#pragma once

#include <vector>
#include <string>

#include "../math/vec3.h"
#include "tmesh.h"

namespace roy {

    namespace mesh {

        class mean_value_coord
        {
        private:
            mean_value_coord(void);

        public:
            template <typename T>
                static void Calcu(const tmesh<T>& fine_mesh, const tmesh<T>& course_mesh, std::vector<std::vector<T> >& W)
            {
                int n = (int) fine_mesh._v_pos.size();
                int m = (int) course_mesh._v_pos.size();

                W.resize(n);

                for (int i = 0; i < n; i++) {
                    CalcuForVertex(fine_mesh._v_pos[i], course_mesh, W[i]);
                }
            }

            template <typename T>
                static bool LoadFromFile_W(const std::string filename, std::vector<std::vector<T> >& W)
            {
                FILE* fp = fopen(filename.c_str(), "r");

                if (fp == NULL) {
                    return false;
                }

                int n, m;
                fscanf(fp, "%d %d", &n, &m);

                char fmt[1024];
                const char* p = get_type_indicator<T>();
                sprintf(fmt, "%%%s", p);
                W.resize(n);
                for (int i = 0; i < n; i++) {
                    W[i].resize(m);
                    for (int j = 0; j < m; j++) {
                        fscanf(fp, fmt, &W[i][j]);
                    }
                }

                fclose(fp);
                return true;
            }

            template <typename T>
                static bool SaveToFile_W(const std::string filename, const std::vector<std::vector<T> >& W)
            {
                FILE* fp = fopen(filename.c_str(), "w");

                if (fp == NULL) {
                    return false;
                }

                int n = (int) W.size(), m = (int) W[0].size();
                fprintf(fp, "%d %d\n", n, m);

                char fmt[1024];
                const char* p = get_type_indicator<T>();
                sprintf(fmt, " %%.8%s", p);
                for (int i = 0; i < n; i++) {
                    for (int j = 0; j < m; j++) {
                        fprintf(fp, fmt, W[i][j]);
                    }
                    fprintf(fp, "\n");
                }

                fclose(fp);
                return true;
            }


            template <typename T>
                static bool LoadFromFile_BinW(const std::string filename, std::vector<std::vector<T> >& W)
            {
                FILE* fp = fopen(filename.c_str(), "rb");

                if (fp == NULL) {
                    return false;
                }

                int n, m;
                fread(&n, sizeof(int), 1, fp);
                fread(&m, sizeof(int), 1, fp);

                W.resize(n);
                for (int i = 0; i < n; i++) {
                    W[i].resize(m);
                    fread(&W[i][0], sizeof(T), m, fp);
                }

                fclose(fp);
                return true;
            }

            template <typename T>
                static bool SaveToFile_BinW(const std::string filename, const std::vector<std::vector<T> >& W)
            {
                FILE* fp = fopen(filename.c_str(), "wb");

                if (fp == NULL) {
                    return false;
                }

                int n = (int) W.size(), m = (int) W[0].size();
                fwrite(&n, sizeof(int), 1, fp);
                fwrite(&m, sizeof(int), 1, fp);

                for (int i = 0; i < n; i++) {
                    assert(W[i].size() == m);
                    fwrite(&W[i][0], sizeof(T), m, fp);
                }

                fclose(fp);
                return true;
            }

            template <typename T>
                static void Reconstructmesh(const tmesh<T>& course_mesh, const std::vector<std::vector<T> >& W, tmesh<T>& recon_mesh)
            {
                int n = (int) W.size();
                int m = (int) course_mesh._v_pos.size();
                recon_mesh._v_cnt = n;
                recon_mesh._v_pos.resize(n);

                for (int i = 0; i < n; i++) {
                    recon_mesh._v_pos[i] = vec3<T>::math::zero();
                    for (int j = 0; j < m; j++) {
                        recon_mesh._v_pos[i] += course_mesh._v_pos[j] * W[i][j];
                    }
                }
            }

            template <typename T>
                static T Error(const tmesh<T>& fine_mesh, const tmesh<T>& course_mesh, const std::vector<std::vector<T> >& W)
            {
                tmesh<T> recon_mesh;
                Reconstructmesh(course_mesh, W, recon_mesh);
                return fine_mesh.Error(recon_mesh);
            }

            template <typename T>
                static T RelativeError(const tmesh<T>& fine_mesh, const tmesh<T>& course_mesh, const std::vector<std::vector<T> >& W)
            {
                tmesh<T> recon_mesh;
                Reconstructmesh(course_mesh, W, recon_mesh);
                return fine_mesh.RelativeError(recon_mesh);
            }

/*
            struct _queue {
                T w;
                int row;
                int col;
                bool operator < (const _queue& another) {
                    return w < another.w;
                }
            };

            void mean_value_coord::RefineWeight(const tmesh& fine_mesh, const tmesh& course_mesh, AAfloat& W, const Afloat& bounds, Afloat& errs, T err_threshold)
            {
                int n = (int) W.size();
                int m = (int) W[0].size();
                int i;

                std::vector<_queue> q(n * m);

                for (i = 0; i < n; i++) {
                    for (int j = 0; j < m; j++) {
                        q[i * m + j].w = fabs(W[i][j]);
                        q[i * m + j].row = i;
                        q[i * m + j].col = j;
                    }
                }

                std::sort(q.begin(), q.end());

                AAfloat t = W;
                int r, c, z_cnt;
                T err;
                int slot = 0;
                errs = Afloat(bounds.size(), 1.0f);

                for (i = 0; i < (int) q.size(); i++) {
                    if (lt(bounds[slot], q[i].w)) {
                        err = Error(fine_mesh, course_mesh, t) / fine_mesh.bounding_sphere_radius;
                        errs[slot++] = err;
                        printf("%.3f deleted -> %.3f%%\n", (T) i / q.size(), err * 100);
                        if (err > err_threshold) {
                            break;
                        }
                    }
                    r = q[i].row;
                    c = q[i].col;

                    z_cnt = 0;
                    for (int j = 0; j < m; j++) {
                        if (math::zero(t[r][j])) {
                            z_cnt++;
                        }
                    }
                    if (z_cnt == m - 1){
                        break;
                    }

                    t[r][c] = 0.0f;
                }
                err = RelativeError(fine_mesh, course_mesh, t);
                errs[slot++] = err;
                printf("%.3f deleted -> %.3f%%\n", (T) (i + 1) / q.size(), err * 100);
                W = t;
            }

            void mean_value_coord::Distribution(const AAfloat& W, const Afloat& bounds, Aint& cnt)
            {
                int slots = (int) bounds.size() - 1;
                cnt = Aint(slots, 0);

                for (int i = 0; i < (int) W.size(); i++) {
                    for (int j = 0; j < (int) W[i].size(); j++) {
                        for (int k = 0; k < slots; k++) {
                            if (!lt(fabs(W[i][j]), bounds[k]) && lt(fabs(W[i][j]), bounds[k + 1])) {
                                cnt[k]++;
                                break;
                            }
                        }
                    }
                }
            }
*/

        //protected:
        public:
            template <typename T>
                static void CalcuForVertex(const math::vec3<T>& x, const tmesh<T>& course_mesh, std::vector<T>& w)
            {
                int m = (int) course_mesh._v_pos.size();

                const std::vector<math::vec3<T> >& p = course_mesh._v_pos;
                const std::vector<std::vector<int> >& face = course_mesh._f_2_v;

                w = std::vector<T>(m);
                std::vector<T> d(m);
                std::vector<math::vec3<T> > u(m);

                math::vec3<T> l;
                math::vec3<T> theta;
                T h;
                math::vec3<T> c;
                math::vec3<T> s;
                T det;
                T sign;


                for (int i = 0; i < m; i++) {
                    d[i] = (T) (p[i] - x);
                    if (math::is_zero(d[i])) {
                        w[i] = 1.0f;
                        return;
                    }
                    u[i] = (p[i] - x).norm_v();
                }

                T w0, w1, w2;
                T total_w = math::zero<T>();

                for (int i = 0; i < (int) course_mesh._f_2_v.size(); i++) {
                    const std::vector<int>& f = face[i];

                    l[0] = (T) (u[f[1]] - u[f[2]]);
                    l[1] = (T) (u[f[2]] - u[f[0]]);
                    l[2] = (T) (u[f[0]] - u[f[1]]);
                    //                std::cout << "l " << l[0] << " " << l[1] << " " << l[2] << std::endl;

                    theta[0] = asin(math::clamp(l[0] / 2, -math::one<T>(), math::one<T>())) * 2;
                    theta[1] = asin(math::clamp(l[1] / 2, -math::one<T>(), math::one<T>())) * 2;
                    theta[2] = asin(math::clamp(l[2] / 2, -math::one<T>(), math::one<T>())) * 2;
                    //                std::cout << "theta " << theta[0] << " " << theta[1] << " " << theta[2] << std::endl;

                    h = (theta[0] + theta[1] + theta[2]) / 2;

                    if (math::is_zero(h - math::pi<T>())) {
                        w0 = sin(theta[0]) * d[f[2]] * d[f[1]];
                        w1 = sin(theta[1]) * d[f[0]] * d[f[2]];
                        w2 = sin(theta[2]) * d[f[1]] * d[f[0]];
                        total_w = w0 + w1 + w2;
                        w = std::vector<T>(m);
                        w[f[0]] = w0 / total_w;
                        w[f[1]] = w1 / total_w;
                        w[f[2]] = w2 / total_w;
                        //                    std::cout << w0 << " " << w1 << " " << w2 << " " << total_w << std::endl;
                        //                    std::cout << w[f[0]] << " " << w[f[0]] << " " << w[f[0]] << std::endl;
                        return;
                    }

                    c[0] = sin(h) * sin(h - theta[0]) * 2 / (sin(theta[1]) * sin(theta[2])) - math::one<T>();
                    c[1] = sin(h) * sin(h - theta[1]) * 2 / (sin(theta[2]) * sin(theta[0])) - math::one<T>();
                    c[2] = sin(h) * sin(h - theta[2]) * 2 / (sin(theta[0]) * sin(theta[1])) - math::one<T>();
                    //                std::cout << "c " << c[0] << " " << c[1] << " " << c[2] << std::endl;

                    det = (u[f[0]] ^ u[f[1]]) * u[f[2]];

                    if (math::is_zero(det)) {
                        continue;
                    }

                    if (math::is_positive(det))
                        sign = math::one<T>();
                    else if (math::is_negative(det))
                        sign = -math::one<T>();
                    else
                        sign = math::zero<T>();

                    s[0] = sign * sqrt(math::clamp(math::one<T>() - math::sq(c[0]), math::zero<T>(), math::one<T>()));
                    s[1] = sign * sqrt(math::clamp(math::one<T>() - math::sq(c[1]), math::zero<T>(), math::one<T>()));
                    s[2] = sign * sqrt(math::clamp(math::one<T>() - math::sq(c[2]), math::zero<T>(), math::one<T>()));
                    //                std::cout << "s " << s[0] << " " << s[1] << " " << s[2] << std::endl;

                    if (math::is_zero(s[0]) || math::is_zero(s[1]) || math::is_zero(s[2])) {
                        continue;
                    }

                    if (math::is_zero(sin(theta[0])) || math::is_zero(sin(theta[1])) || math::is_zero(sin(theta[2]))) {
                        continue;
                    }

                    w0 = (theta[0] - c[1] * theta[2] - c[2] * theta[1]) / (d[f[0]] * sin(theta[1]) * s[2]);
                    w1 = (theta[1] - c[2] * theta[0] - c[0] * theta[2]) / (d[f[1]] * sin(theta[2]) * s[0]);
                    w2 = (theta[2] - c[0] * theta[1] - c[1] * theta[0]) / (d[f[2]] * sin(theta[0]) * s[1]);
                    w[f[0]] += w0;
                    w[f[1]] += w1;
                    w[f[2]] += w2;
                    //                std::cout << "w " << w0 << " " << w1 << " " << w2 << std::endl;
                }

                // Unify the weights.
                total_w = math::zero<T>();
                for (int i = 0; i < m; i++) {
                    total_w += w[i];
                }
                for (int i = 0; i < m; i++) {
                    w[i] /= total_w;
                }
            }

        };

    }

}
